############################################################
############         Appendix Figure A2         ############
############################################################

# Code is in "SummaryStats_and_Balance.R"


############################################################
############         Appendix Figure A3         ############
############################################################

# Code is in "SummaryStats_and_Balance.R"


############################################################
############         Appendix Figure A4         ############
############################################################

hypertension <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #hypertension_9; DROP TABLE IF EXISTS #hypertension_10;
SELECT DISTINCT ICD9SID INTO #hypertension_9 FROM OMHSP_PERC_Share.DOEx.Lookup_ICD9 WHERE EH_HYPERTENS=1;
SELECT DISTINCT ICD10SID INTO #hypertension_10 FROM OMHSP_PERC_Share.DOEx.Lookup_ICD10 WHERE EH_HYPERTENS=1

SELECT cohort.PatientICN
FROM CDWWork.Outpat.VDiagnosis AS diag
INNER JOIN #hypertension_9
  ON diag.ICD9SID=#hypertension_9.ICD9SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE diag.VisitDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2014-10-01 00:00:00' AS datetime2(0))
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Outpat.VDiagnosis AS diag
INNER JOIN #hypertension_10
  ON diag.ICD10SID=#hypertension_10.ICD10SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE diag.VisitDateTime>=CAST('2014-10-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Inpat.InpatientDiagnosis AS diag
INNER JOIN #hypertension_9
  ON diag.ICD9SID=#hypertension_9.ICD9SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND DischargeDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE DischargeDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND DischargeDateTime<CAST('2014-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Inpat.InpatientDiagnosis AS diag
INNER JOIN #hypertension_10
  ON diag.ICD10SID=#hypertension_10.ICD10SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND DischargeDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE DischargeDateTime>=CAST('2014-10-01 00:00:00' AS datetime2(0)) AND DischargeDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
)
")

diabetes <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; 
DROP TABLE IF EXISTS #diabetes_9; DROP TABLE IF EXISTS #diabetes_10;
SELECT DISTINCT ICD9SID INTO #diabetes_9 FROM OMHSP_PERC_Share.DOEx.Lookup_ICD9 WHERE EH_UNCDIAB=1 OR EH_COMDIAB=1
SELECT DISTINCT ICD10SID INTO #diabetes_10 FROM OMHSP_PERC_Share.DOEx.Lookup_ICD10 WHERE EH_UNCDIAB=1 OR EH_COMDIAB=1

SELECT cohort.PatientICN
FROM CDWWork.Outpat.VDiagnosis AS diag
INNER JOIN #diabetes_9
  ON diag.ICD9SID=#diabetes_9.ICD9SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE diag.VisitDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2014-10-01 00:00:00' AS datetime2(0))
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Outpat.VDiagnosis AS diag
INNER JOIN #diabetes_10
  ON diag.ICD10SID=#diabetes_10.ICD10SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE diag.VisitDateTime>=CAST('2014-10-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Inpat.InpatientDiagnosis AS diag
INNER JOIN #diabetes_9
  ON diag.ICD9SID=#diabetes_9.ICD9SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND DischargeDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE DischargeDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND DischargeDateTime<CAST('2014-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Inpat.InpatientDiagnosis AS diag
INNER JOIN #diabetes_10
  ON diag.ICD10SID=#diabetes_10.ICD10SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND DischargeDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE DischargeDateTime>=CAST('2014-10-01 00:00:00' AS datetime2(0)) AND DischargeDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
)
")

cholesterol <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; 
DROP TABLE IF EXISTS #cholesterol_9; DROP TABLE IF EXISTS #cholesterol_10;
SELECT DISTINCT ICD9SID INTO #cholesterol_9 FROM OMHSP_PERC_Share.DOEx.Lookup_ICD9 WHERE SUBSTRING(ICD9Code, 1, 5) IN ('272.0', '272.1', '272.2', '272.5') 
SELECT DISTINCT ICD10SID INTO #cholesterol_10 FROM OMHSP_PERC_Share.DOEx.Lookup_ICD10 WHERE SUBSTRING(ICD10Code, 1, 5) IN ('E78.0', 'E78.1', 'E78.2', 'E78.6')

SELECT cohort.PatientICN
FROM CDWWork.Outpat.VDiagnosis AS diag
INNER JOIN #cholesterol_9
  ON diag.ICD9SID=#cholesterol_9.ICD9SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE diag.VisitDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2014-10-01 00:00:00' AS datetime2(0))
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Outpat.VDiagnosis AS diag
INNER JOIN #cholesterol_10
  ON diag.ICD10SID=#cholesterol_10.ICD10SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE diag.VisitDateTime>=CAST('2014-10-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Inpat.InpatientDiagnosis AS diag
INNER JOIN #cholesterol_9
  ON diag.ICD9SID=#cholesterol_9.ICD9SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND DischargeDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE DischargeDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND DischargeDateTime<CAST('2014-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT cohort.PatientICN
FROM CDWWork.Inpat.InpatientDiagnosis AS diag
INNER JOIN #cholesterol_10
  ON diag.ICD10SID=#cholesterol_10.ICD10SID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON diag.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND DischargeDateTime BETWEEN DATEADD(day, -365, AppointmentRequestDate) AND DATEADD(day, 30, AppointmentRequestDate)
WHERE DischargeDateTime>=CAST('2014-10-01 00:00:00' AS datetime2(0)) AND DischargeDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
)
")

pcp_claims$PriorHypertension <- (pcp_claims$PatientICN %in% hypertension$PatientICN)
pcp_claims$PriorDiabetes <- (pcp_claims$PatientICN %in% diabetes$PatientICN)
pcp_claims$PriorCholesterol <- (pcp_claims$PatientICN %in% cholesterol$PatientICN)


a1c <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #a1c_lab; DROP TABLE IF EXISTS #a1c_loinc
SELECT DISTINCT a.LabChemTestSID INTO #a1c_lab 
FROM CDWWork.Dim.LabChemTest AS a
INNER JOIN CDWWork.Dim.NationalVALabCode AS b
  ON a.NationalVALabCodeSID=b.NationalVALabCodeSID
WHERE LabProcedure LIKE '%a1c%' OR LabChemPrintTestName LIKE '%a1c%' OR LabChemTestName LIKE '%a1c%'

SELECT DISTINCT LOINCSID INTO #a1c_loinc 
FROM CDWWork.Dim.LOINC
WHERE LOINC IN ('4548-4', '17855-8', '17856-6', '59261-8', '86910-7', '41995-2', '4549-2', '62388-4', '71875-9')

SELECT
    cohort.PatientICN, CAST(chem.LabChemSpecimenDateTime AS DATE) AS a1cDate, LabChemResultNumericValue AS a1c
  FROM CDWWork.Chem.PatientLabChem AS chem
  INNER JOIN #a1c_lab
    ON chem.LabChemTestSID=#a1c_lab.LabChemTestSID
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON chem.PatientSID=pat.PatientSID
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON pat.PatientICN=cohort.PatientICN AND LabChemSpecimenDateTime BETWEEN DATEADD(day, -360, AppointmentRequestDate) AND DATEADD(day, 365*3, AppointmentRequestDate)
  WHERE chem.LabChemSpecimenDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
  AND LabChemResultNumericValue<20
  UNION(
  SELECT
    cohort.PatientICN, CAST(chem.LabChemSpecimenDateTime AS DATE) AS a1cDate, LabChemResultNumericValue AS a1c
  FROM CDWWork.Chem.PatientLabChem AS chem
  INNER JOIN #a1c_loinc
    ON chem.LOINCSID=#a1c_loinc.LOINCSID
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON chem.PatientSID=pat.PatientSID
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON pat.PatientICN=cohort.PatientICN AND LabChemSpecimenDateTime BETWEEN DATEADD(day, -360, AppointmentRequestDate) AND DATEADD(day, 365*3, AppointmentRequestDate)
  WHERE chem.LabChemSpecimenDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
  AND LabChemResultNumericValue<20
) 
")

a1c <- data.table(a1c)
a1c$a1cDate <- as.Date(a1c$a1cDate)
a1c <- a1c[, .(a1c=median(a1c)), by=c("PatientICN", "a1cDate")]
a1c <- merge(x=a1c, y=pcp_claims[,c("PatientICN", "AppointmentRequestDate")], by.x="PatientICN", by.y="PatientICN")
a1c$AppointmentRequestDate <- as.Date(a1c$AppointmentRequestDate)
a1c$RelativeHalf <- floor(as.numeric(a1c$a1cDate-a1c$AppointmentRequestDate)/180)

a1c <- a1c[RelativeHalf<=5, .(a1c_median=median(a1c), a1c_mean=mean(a1c)), by=c("PatientICN", "RelativeHalf")]
a1c <- merge(x=a1c, y=pcp_claims[,c("PatientICN", "PriorDiabetes", "DeathRelativeMonth", "FE_MH_std", "FE_PH_std", "FE_ACSC_std", "RaceEthnicity", "Married", "Medicare", "Medicaid", "PeriodOfService",
                                    "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount", "PriorVisit", "MH_dummy_lag1_v1", 
                                    "PH_dummy_lag1_v1", "ACSC_lag1", "WaitTime", "Age_bins", "EnrollmentPriority", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")], by.x="PatientICN", by.y="PatientICN")

a1c$RelativeHalf <- as.factor(a1c$RelativeHalf); a1c$RelativeHalf <- relevel(a1c$RelativeHalf, "-1")

a1c$a1cControl <- (a1c$a1c_median<5.7)


BP <-  sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #pc_stopcode;
SELECT DISTINCT StopCodeSID INTO #pc_stopcode FROM CDWWork.Dim.StopCode WHERE StopCode IN (301, 322, 323, 348, 350, 342)

SELECT
    cohort.PatientICN, CAST(VitalSignTakenDateTime AS DATE) AS VitalDate, Systolic, Diastolic
  FROM CDWWork.Vital.VitalSign AS pain
  INNER JOIN CDWWork.Dim.Location AS loc
    ON pain.LocationSID=loc.LocationSID
  LEFT OUTER JOIN #pc_stopcode AS psc
    ON loc.PrimaryStopCodeSID=psc.StopCodeSID
  LEFT OUTER JOIN #pc_stopcode AS ssc
    ON loc.SecondaryStopCodeSID=ssc.StopCodeSID
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON pain.PatientSID=pat.PatientSID
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON pat.PatientICN=cohort.PatientICN AND VitalSignTakenDateTime BETWEEN DATEADD(day, -360, AppointmentRequestDate) AND DATEADD(day, 365*3, AppointmentRequestDate)
  WHERE Systolic>=50 AND Systolic<=250 AND Diastolic>=50 AND Diastolic<=250
  AND (psc.StopCodeSID IS NOT NULL OR ssc.StopCodeSID IS NOT NULL)
  AND VitalSignTakenDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND VitalSignTakenDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
")
BP <- data.table(BP)

BP$VitalDate <- as.Date(BP$VitalDate)
BP <- BP[, .(Systolic=median(Systolic), Diastolic=median(Diastolic)), by=c("PatientICN", "VitalDate")]
BP <- merge(x=BP, y=pcp_claims[,c("PatientICN", "AppointmentRequestDate")], by.x="PatientICN", by.y="PatientICN")
BP$AppointmentRequestDate <- as.Date(BP$AppointmentRequestDate)
BP$RelativeHalf <- floor(as.numeric(BP$VitalDate-BP$AppointmentRequestDate)/180)

BP <- BP[RelativeHalf<=5, .(Systolic=median(Systolic), Diastolic=median(Diastolic)), by=c("PatientICN", "RelativeHalf")]
BP <- merge(x=BP, y=pcp_claims[,c("PatientICN", "PriorHypertension", "PriorDiabetes", "DeathRelativeMonth", "FE_MH_std", "FE_PH_std", "FE_ACSC_std", "RaceEthnicity", "Married", "Medicare", "Medicaid", "PeriodOfService",
                                  "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount", "PriorVisit", "MH_dummy_lag1_v1", 
                                  "PH_dummy_lag1_v1", "ACSC_lag1", "WaitTime", "Age_bins", "EnrollmentPriority", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")], by.x="PatientICN", by.y="PatientICN")

BP$RelativeHalf <- as.factor(BP$RelativeHalf); BP$RelativeHalf <- relevel(BP$RelativeHalf, "-1")

BP$SystolicControl <- (BP$Systolic<140); BP$DiastolicControl <- (BP$Diastolic<90)


# MPR among hypertensive patients
MPR <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #hypertensives;

SELECT DISTINCT DrugNameWithoutDose INTO #hypertensives
FROM CDWWork.Dim.DrugNameWithoutDose AS without
INNER JOIN CDWWork.Dim.NationalDrug AS nat
  ON without.DrugNameWithoutDoseSID=nat.DrugNameWithoutDoseSID
INNER JOIN CDWWork.Dim.DrugClass AS class
  ON nat.PrimaryDrugClassSID=class.DrugClassSID
WHERE DrugClassCode IN ('CV800', 'CV805', 'CV806', 'CV100', 'CV700', 'CV200');

SELECT
  cohort.PatientICN, CAST(TrialStartDateTime AS date) AS TrialStartDate, MPR_Trial AS MPR
  FROM OIT_ROCKIES.DOEx.OIT_Rockies_MPR_Drug AS mpr
  INNER JOIN #hypertensives
    ON mpr.DrugNameWithoutDose=#hypertensives.DrugNameWithoutDose
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON mpr.PatientICN=cohort.PatientICN AND TrialStartDateTime BETWEEN DATEADD(day, -360, AppointmentRequestDate) AND DATEADD(day, 365*3, AppointmentRequestDate)
  WHERE TrialStartDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND TrialStartDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
AND MonthsInTreatment>0
")
MPR <- data.table(MPR)

MPR$TrialStartDate <- as.Date(MPR$TrialStartDate)
MPR$MPR[which(MPR$MPR>1)] <- 1
MPR <- merge(x=MPR, y=pcp_claims[,c("PatientICN", "AppointmentRequestDate")], by.x="PatientICN", by.y="PatientICN")
MPR$AppointmentRequestDate <- as.Date(MPR$AppointmentRequestDate)
MPR$RelativeHalf <- floor(as.numeric(MPR$TrialStartDate-MPR$AppointmentRequestDate)/180)


MPR <- MPR[RelativeHalf<=5, .(MPR_mean=mean(MPR)), by=c("PatientICN", "RelativeHalf")]
MPR <- merge(x=MPR, y=pcp_claims[,c("PatientICN", "PriorHypertension", "PriorDiabetes", "DeathRelativeMonth", "FE_MH_std", "FE_PH_std", "FE_ACSC_std", "RaceEthnicity", "Married", "Medicare", "Medicaid", "PeriodOfService",
                                    "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount", "PriorVisit", "MH_dummy_lag1_v1", 
                                    "PH_dummy_lag1_v1", "ACSC_lag1", "WaitTime", "Age_bins", "EnrollmentPriority", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")], by.x="PatientICN", by.y="PatientICN")

MPR$RelativeHalf <- as.factor(MPR$RelativeHalf)
MPR$RelativeHalf <- relevel(MPR$RelativeHalf, "-1")

MPR$MPRHigh <- (MPR$MPR_mean>=0.8)

ldl <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #ldl_lab; DROP TABLE IF EXISTS #ldl_loinc
SELECT DISTINCT a.LabChemTestSID INTO #ldl_lab 
FROM CDWWork.Dim.LabChemTest AS a
INNER JOIN CDWWork.Dim.NationalVALabCode AS b
  ON a.NationalVALabCodeSID=b.NationalVALabCodeSID
WHERE LabProcedure LIKE '%LDL%' OR LabChemPrintTestName LIKE '%LDL%' OR LabChemTestName LIKE '%LDL%'

SELECT DISTINCT LOINCSID INTO #ldl_loinc 
FROM CDWWork.Dim.LOINC
WHERE Component LIKE '%CHOLESTEROL.IN LDL%'

SELECT
    cohort.PatientICN, CAST(chem.LabChemSpecimenDateTime AS DATE) AS ldlDate, LabChemResultNumericValue AS LDL, Units
  FROM CDWWork.Chem.PatientLabChem AS chem
  INNER JOIN #ldl_lab
    ON chem.LabChemTestSID=#ldl_lab.LabChemTestSID
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON chem.PatientSID=pat.PatientSID
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON pat.PatientICN=cohort.PatientICN AND LabChemSpecimenDateTime BETWEEN DATEADD(day, -360, AppointmentRequestDate) AND DATEADD(day, 365*3, AppointmentRequestDate)
  WHERE chem.LabChemSpecimenDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0)) 
  AND Units IN ('mg/dl', 'mg/dL', 'mG/dL', 'MG/DL', 'mg/dl')
  UNION(
  SELECT
    cohort.PatientICN, CAST(chem.LabChemSpecimenDateTime AS DATE) AS ldlDate, LabChemResultNumericValue AS LDL, Units
  FROM CDWWork.Chem.PatientLabChem AS chem
  INNER JOIN #ldl_loinc
    ON chem.LOINCSID=#ldl_loinc.LOINCSID
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON chem.PatientSID=pat.PatientSID
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON pat.PatientICN=cohort.PatientICN AND LabChemSpecimenDateTime BETWEEN DATEADD(day, -360, AppointmentRequestDate) AND DATEADD(day, 365*3, AppointmentRequestDate)
  WHERE chem.LabChemSpecimenDateTime>=CAST('2003-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2021-01-01 00:00:00' AS datetime2(0))
  AND Units IN ('mg/dl', 'mg/dL', 'mG/dL', 'MG/DL', 'mg/dl')
) 
")

ldl <- data.table(ldl)
ldl$ldlDate <- as.Date(ldl$ldlDate)
ldl <- ldl[LDL>=0 & LDL<=400]

ldl <- ldl[, .(LDL=median(LDL)), by=c("PatientICN", "ldlDate")]
ldl <- merge(x=ldl, y=pcp_claims[,c("PatientICN", "AppointmentRequestDate")], by.x="PatientICN", by.y="PatientICN")
ldl$AppointmentRequestDate <- as.Date(ldl$AppointmentRequestDate)
ldl$RelativeHalf <- floor(as.numeric(ldl$ldlDate-ldl$AppointmentRequestDate)/180)

ldl <- ldl[RelativeHalf<=5, .(LDL_mean=mean(LDL)), by=c("PatientICN", "RelativeHalf")]
ldl <- merge(x=ldl, y=pcp_claims[,c("PatientICN", "PriorHypertension", "PriorDiabetes", "PriorCholesterol", "DeathRelativeMonth", "FE_MH_std", "FE_PH_std", "FE_ACSC_std", "RaceEthnicity", "Married", "Medicare", "Medicaid", "PeriodOfService",
                                    "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount", "PriorVisit", "MH_dummy_lag1_v1", 
                                    "PH_dummy_lag1_v1", "ACSC_lag1", "WaitTime", "Age_bins", "EnrollmentPriority", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")], by.x="PatientICN", by.y="PatientICN")

ldl$RelativeHalf <- as.factor(ldl$RelativeHalf); ldl$RelativeHalf <- relevel(ldl$RelativeHalf, "-1")

ldl$LDLControl <- (ldl$LDL_mean<160)


##### PLOT (save as 18x20)
par(mfrow=c(4,3), mar=c(4.5,4.7,3,2.5), oma=c(1,3,1,1))
### A1c
dx_control <- felm(a1cControl~FE_MH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=a1c[PriorDiabetes==1 & DeathRelativeMonth>36])
dx_control <- rbind(summary(dx_control)$coef[35, 1:2], c(0, NA), summary(dx_control)$coef[36:41, 1:2])
plot(y=dx_control[,1], x=-2:5, ylim=c(-0.01,0.015), ylab="", xlab="Half-Year Relative to Assignment", main="Mental Health Effectiveness", cex.main=2.5, cex.axis=1.5, cex.lab=1.5, font.main=1, pch=19, cex=1.5)
mtext(text="HbA1c Control", 2, line=5, cex=1.8)
arrows(x0=-2:5, x1=-2:5, y0=dx_control[,1]-1.96*dx_control[,2], y1=dx_control[,1]+1.96*dx_control[,2], length=0, angle=90, lwd=2); abline(h=0)

dx_control <- felm(a1cControl~FE_PH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=a1c[PriorDiabetes==1 & DeathRelativeMonth>36])
dx_control <- rbind(summary(dx_control)$coef[35, 1:2], c(0, NA), summary(dx_control)$coef[36:41, 1:2])
plot(y=dx_control[,1], x=-2:5, ylim=c(-0.01,0.015), ylab="", xlab="Half-Year Relative to Assignment", main="Circulatory Condition Effectiveness", cex.main=2.5, cex.axis=1.5, cex.lab=1.5, font.main=1, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=dx_control[,1]-1.96*dx_control[,2], y1=dx_control[,1]+1.96*dx_control[,2], length=0, angle=90, lwd=2); abline(h=0)

dx_control <- felm(a1cControl~FE_ACSC_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=a1c[PriorDiabetes==1 & DeathRelativeMonth>36])
dx_control <- rbind(summary(dx_control)$coef[35, 1:2], c(0, NA), summary(dx_control)$coef[36:41, 1:2])
plot(y=dx_control[,1], x=-2:5, ylim=c(-0.01,0.015), ylab="", xlab="Half-Year Relative to Assignment", main="ACSC Effectiveness", cex.main=2.5, cex.axis=1.5, cex.lab=1.5, font.main=1, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=dx_control[,1]-1.96*dx_control[,2], y1=dx_control[,1]+1.96*dx_control[,2], length=0, angle=90, lwd=2); abline(h=0)


### BP
bp_control <- felm(SystolicControl~FE_MH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=BP[PriorHypertension==1 & DeathRelativeMonth>36])
bp_control <- rbind(summary(bp_control)$coef[35, 1:2], c(0, NA), summary(bp_control)$coef[36:41, 1:2])
plot(y=bp_control[,1], x=-2:5, ylim=c(-0.012,0.015), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
mtext(text="Systolic BP Control", 2, line=5, cex=1.8)
arrows(x0=-2:5, x1=-2:5, y0=bp_control[,1]-1.96*bp_control[,2], y1=bp_control[,1]+1.96*bp_control[,2], length=0, angle=90, lwd=2); abline(h=0)

bp_control <- felm(SystolicControl~FE_PH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=BP[PriorHypertension==1 & DeathRelativeMonth>36])
bp_control <- rbind(summary(bp_control)$coef[35, 1:2], c(0, NA), summary(bp_control)$coef[36:41, 1:2])
plot(y=bp_control[,1], x=-2:5, ylim=c(-0.012,0.015), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=bp_control[,1]-1.96*bp_control[,2], y1=bp_control[,1]+1.96*bp_control[,2], length=0, angle=90, lwd=2); abline(h=0)

bp_control <- felm(SystolicControl~FE_ACSC_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=BP[PriorHypertension==1 & DeathRelativeMonth>36])
bp_control <- rbind(summary(bp_control)$coef[35, 1:2], c(0, NA), summary(bp_control)$coef[36:41, 1:2])
plot(y=bp_control[,1], x=-2:5, ylim=c(-0.012,0.015), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=bp_control[,1]-1.96*bp_control[,2], y1=bp_control[,1]+1.96*bp_control[,2], length=0, angle=90, lwd=2); abline(h=0)


### MPR
mprhigh <- felm(MPRHigh~FE_MH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=MPR[PriorHypertension==1 & DeathRelativeMonth>36])
mprhigh <- rbind(summary(mprhigh)$coef[35, 1:2], c(0, NA), summary(mprhigh)$coef[36:41, 1:2])

plot(y=mprhigh[,1], x=-2:5, ylim=c(-0.015,0.03), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
mtext(text="Medication Adherence", 2, line=5, cex=1.8)
arrows(x0=-2:5, x1=-2:5, y0=mprhigh[,1]-1.96*mprhigh[,2], y1=mprhigh[,1]+1.96*mprhigh[,2], length=0, angle=90, lwd=2); abline(h=0)

mprhigh <- felm(MPRHigh~FE_PH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=MPR[PriorHypertension==1 & DeathRelativeMonth>36])
mprhigh <- rbind(summary(mprhigh)$coef[35, 1:2], c(0, NA), summary(mprhigh)$coef[36:41, 1:2])
plot(y=mprhigh[,1], x=-2:5, ylim=c(-0.015,0.03), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=mprhigh[,1]-1.96*mprhigh[,2], y1=mprhigh[,1]+1.96*mprhigh[,2], length=0, angle=90, lwd=2); abline(h=0)

mprhigh <- felm(MPRHigh~FE_ACSC_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=MPR[PriorHypertension==1 & DeathRelativeMonth>36])
mprhigh <- rbind(summary(mprhigh)$coef[35, 1:2], c(0, NA), summary(mprhigh)$coef[36:41, 1:2])
plot(y=mprhigh[,1], x=-2:5, ylim=c(-0.015,0.03), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=mprhigh[,1]-1.96*mprhigh[,2], y1=mprhigh[,1]+1.96*mprhigh[,2], length=0, angle=90, lwd=2); abline(h=0)


### LDL
ldlmodel <- felm(LDLControl~FE_MH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=ldl[PriorCholesterol==1 & DeathRelativeMonth>36])
ldlmodel <- rbind(summary(ldlmodel)$coef[35, 1:2], c(0, NA), summary(ldlmodel)$coef[36:41, 1:2])
plot(y=ldlmodel[,1], x=-2:5, ylim=c(-0.015,0.02), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
mtext(text="LDL Control", 2, line=5, cex=1.8)
arrows(x0=-2:5, x1=-2:5, y0=ldlmodel[,1]-1.96*ldlmodel[,2], y1=ldlmodel[,1]+1.96*ldlmodel[,2], length=0, angle=90, lwd=2); abline(h=0)

ldlmodel <- felm(LDLControl~FE_PH_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=ldl[PriorCholesterol==1 & DeathRelativeMonth>36])
ldlmodel <- rbind(summary(ldlmodel)$coef[35, 1:2], c(0, NA), summary(ldlmodel)$coef[36:41, 1:2])
plot(y=ldlmodel[,1], x=-2:5, ylim=c(-0.015,0.02), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=ldlmodel[,1]-1.96*ldlmodel[,2], y1=ldlmodel[,1]+1.96*ldlmodel[,2], length=0, angle=90, lwd=2); abline(h=0)

ldlmodel <- felm(LDLControl~FE_ACSC_std*RelativeHalf+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=ldl[PriorCholesterol==1 & DeathRelativeMonth>36])
ldlmodel <- rbind(summary(ldlmodel)$coef[35, 1:2], c(0, NA), summary(ldlmodel)$coef[36:41, 1:2])
plot(y=ldlmodel[,1], x=-2:5, ylim=c(-0.015,0.02), ylab="", xlab="Half-Year Relative to Assignment", main="", cex.main=1.5, cex.axis=1.5, cex.lab=1.5, pch=19, cex=1.5)
arrows(x0=-2:5, x1=-2:5, y0=ldlmodel[,1]-1.96*ldlmodel[,2], y1=ldlmodel[,1]+1.96*ldlmodel[,2], length=0, angle=90, lwd=2); abline(h=0)


############################################################
############         Appendix Table A1         #############
############################################################

### CONSTRUCT EACH CONDITION ON THREE DISJOINT SUBSAMPLES
set.seed(1)
pcp_claims <- pcp_claims[sample(1:(dim(pcp_claims)[1])),]
pcp_claims$ConstructionSample <- rep(c("MH", "PH", "ACSC"), length.out=dim(pcp_claims)[1])

MHsample <- pcp_claims[ConstructionSample=="MH"]
PHsample <- pcp_claims[ConstructionSample=="PH"]
ACSCsample <- pcp_claims[ConstructionSample=="ACSC"]

MHsample$residual_MH <- resid(felm(MH_dummy_cum3_v1~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=MHsample))
PHsample$residual_PH <- resid(felm(PH_dummy_cum3_v1~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=PHsample))
ACSCsample$residual_ACSC <- resid(felm(ACSC_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=ACSCsample))

MHdrift <- MHsample[, `:=` (Sum_Resid_MH=sum(residual_MH, na.rm=T), AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]
PHdrift <- PHsample[, `:=` (Sum_Resid_PH=sum(residual_PH, na.rm=T),  AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]
ACSCdrift <- ACSCsample[, `:=` (Sum_Resid_ACSC=sum(residual_ACSC, na.rm=T), AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]

MHdrift$Z_MH_year <- (MHdrift$Sum_Resid_MH - MHdrift$residual_MH)/(MHdrift$AnnualCaseCount - 1)
PHdrift$Z_PH_year <- (PHdrift$Sum_Resid_PH - PHdrift$residual_PH)/(PHdrift$AnnualCaseCount - 1)
ACSCdrift$Z_ACSC_year <- (ACSCdrift$Sum_Resid_ACSC - ACSCdrift$residual_ACSC)/(ACSCdrift$AnnualCaseCount - 1)

MHdrift <- MHdrift[,c("PatientICN", "MH_dummy_cum3_v1", "Z_MH_year", "StaffSSN", "AnnualCaseCount", "AssignYear", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")]
PHdrift <- PHdrift[,c("PatientICN", "PH_dummy_cum3_v1", "Z_PH_year", "StaffSSN", "AnnualCaseCount", "AssignYear", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")]
ACSCdrift <- ACSCdrift[,c("PatientICN", "ACSC_cum3", "Z_ACSC_year", "StaffSSN", "AnnualCaseCount", "AssignYear", "InstitutionSID", "YearXMonth", "DayOfWeek", "TimeOut_bins")]

MHdrift$N_bin <- "1"; MHdrift$N_bin[which(MHdrift$AnnualCaseCount>=10 & MHdrift$AnnualCaseCount<25)] <- "2"; MHdrift$N_bin[which(MHdrift$AnnualCaseCount>=25 & MHdrift$AnnualCaseCount<50)] <- "3"; MHdrift$N_bin[which(MHdrift$AnnualCaseCount>=50)] <- "4"
PHdrift$N_bin <- "1"; PHdrift$N_bin[which(PHdrift$AnnualCaseCount>=10 & PHdrift$AnnualCaseCount<25)] <- "2"; PHdrift$N_bin[which(PHdrift$AnnualCaseCount>=25 & PHdrift$AnnualCaseCount<50)] <- "3"; PHdrift$N_bin[which(PHdrift$AnnualCaseCount>=50)] <- "4"
ACSCdrift$N_bin <- "1"; ACSCdrift$N_bin[which(ACSCdrift$AnnualCaseCount>=10 & ACSCdrift$AnnualCaseCount<25)] <- "2"; ACSCdrift$N_bin[which(ACSCdrift$AnnualCaseCount>=25 & ACSCdrift$AnnualCaseCount<50)] <- "3"; ACSCdrift$N_bin[which(ACSCdrift$AnnualCaseCount>=50)] <- "4"

MHdrift_wide <- reshape(MHdrift, v.names=c("Z_MH_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide"); MHdrift_wide <- data.frame(MHdrift_wide); MHdrift_wide[is.na(MHdrift_wide)] <- 0; MHdrift_wide <- data.table(MHdrift_wide)
PHdrift_wide <- reshape(PHdrift, v.names=c("Z_PH_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide"); PHdrift_wide <- data.frame(PHdrift_wide); PHdrift_wide[is.na(PHdrift_wide)] <- 0; PHdrift_wide <- data.table(PHdrift_wide)
ACSCdrift_wide <- reshape(ACSCdrift, v.names=c("Z_ACSC_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide"); ACSCdrift_wide <- data.frame(ACSCdrift_wide); ACSCdrift_wide[is.na(ACSCdrift_wide)] <- 0; ACSCdrift_wide <- data.table(ACSCdrift_wide)

MHdrift <- merge(x=MHdrift, y=MHdrift_wide[,-c("PatientICN", "MH_dummy_cum3_v1")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)
PHdrift <- merge(x=PHdrift, y=PHdrift_wide[,-c("PatientICN", "PH_dummy_cum3_v1")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)
ACSCdrift <- merge(x=ACSCdrift, y=ACSCdrift_wide[,-c("PatientICN", "ACSC_cum3")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

MHdrift$FE_MH <- predict(feols(MH_dummy_cum3_v1~N_bin.2005*Z_MH_year.2005+N_bin.2006*Z_MH_year.2006+N_bin.2007*Z_MH_year.2007+N_bin.2008*Z_MH_year.2008+N_bin.2009*Z_MH_year.2009+N_bin.2010*Z_MH_year.2010+N_bin.2011*Z_MH_year.2011+N_bin.2012*Z_MH_year.2012+N_bin.2013*Z_MH_year.2013+N_bin.2014*Z_MH_year.2014+N_bin.2015*Z_MH_year.2015+N_bin.2016*Z_MH_year.2016+N_bin.2017*Z_MH_year.2017|InstitutionSID.x+YearXMonth.x+DayOfWeek.x+TimeOut_bins.x, data=MHdrift))
PHdrift$FE_PH <- predict(feols(PH_dummy_cum3_v1~N_bin.2005*Z_PH_year.2005+N_bin.2006*Z_PH_year.2006+N_bin.2007*Z_PH_year.2007+N_bin.2008*Z_PH_year.2008+N_bin.2009*Z_PH_year.2009+N_bin.2010*Z_PH_year.2010+N_bin.2011*Z_PH_year.2011+N_bin.2012*Z_PH_year.2012+N_bin.2013*Z_PH_year.2013+N_bin.2014*Z_PH_year.2014+N_bin.2015*Z_PH_year.2015+N_bin.2016*Z_PH_year.2016+N_bin.2017*Z_PH_year.2017|InstitutionSID.x+YearXMonth.x+DayOfWeek.x+TimeOut_bins.x, data=PHdrift))
ACSCdrift$FE_ACSC <- predict(feols(ACSC_cum3~N_bin.2005*Z_ACSC_year.2005+N_bin.2006*Z_ACSC_year.2006+N_bin.2007*Z_ACSC_year.2007+N_bin.2008*Z_ACSC_year.2008+N_bin.2009*Z_ACSC_year.2009+N_bin.2010*Z_ACSC_year.2010+N_bin.2011*Z_ACSC_year.2011+N_bin.2012*Z_ACSC_year.2012+N_bin.2013*Z_ACSC_year.2013+N_bin.2014*Z_ACSC_year.2014+N_bin.2015*Z_ACSC_year.2015+N_bin.2016*Z_ACSC_year.2016+N_bin.2017*Z_ACSC_year.2017|InstitutionSID.x+YearXMonth.x+DayOfWeek.x+TimeOut_bins.x, data=ACSCdrift))


MHsample <- merge(x=MHsample, y=MHdrift[,c("PatientICN", "FE_MH")], by.x="PatientICN", by.y="PatientICN", all.x=T)
PHsample <- merge(x=PHsample, y=PHdrift[,c("PatientICN", "FE_PH")], by.x="PatientICN", by.y="PatientICN", all.x=T)
ACSCsample <- merge(x=ACSCsample, y=ACSCdrift[,c("PatientICN", "FE_ACSC")], by.x="PatientICN", by.y="PatientICN", all.x=T)

# Standardize, then merge with the left-out sample
MHsample <- MHsample[, .(FE_MH_subsample=mean(FE_MH.y)), by="StaffSSN"]
MHsample$FE_MH_subsample <- -scale(MHsample$FE_MH_subsample)
MHsample <- merge(x=pcp_claims[ConstructionSample %in% c("PH", "ACSC")], y=MHsample, by.x="StaffSSN", by.y="StaffSSN")

mh <- felm(MH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=MHsample)
ph <- felm(PH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=MHsample)
acsc <- felm(ACSC_cum3~FE_MH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=MHsample)

a <- stargazer(mh, ph, acsc,
          apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),mdigits=2,
          add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(mh$fitted.values),4), 100*round(mean(ph$fitted.values),4), 100*round(mean(acsc$fitted.values),4))))

PHsample <- PHsample[, .(FE_PH_subsample=mean(FE_PH.y)), by="StaffSSN"]
PHsample$FE_MH_subsample <- -scale(PHsample$FE_PH_subsample)
PHsample <- merge(x=pcp_claims[ConstructionSample %in% c("MH", "ACSC")], y=PHsample, by.x="StaffSSN", by.y="StaffSSN")

mh <- felm(MH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=PHsample)
ph <- felm(PH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=PHsample)
acsc <- felm(ACSC_cum3~FE_PH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=PHsample)

b <- stargazer(mh, ph, acsc,
          apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),mdigits=2,
          add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(mh$fitted.values),4), 100*round(mean(ph$fitted.values),4), 100*round(mean(acsc$fitted.values),4))))

ACSCsample <- ACSCsample[, .(FE_ACSC_subsample=mean(FE_ACSC.y)), by="StaffSSN"]
ACSCsample$FE_MH_subsample <- -scale(ACSCsample$FE_ACSC_subsample)
ACSCsample <- merge(x=pcp_claims[ConstructionSample %in% c("PH", "MH")], y=ACSCsample, by.x="StaffSSN", by.y="StaffSSN")

mh <- felm(MH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=ACSCsample)
ph <- felm(PH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=ACSCsample)
acsc <- felm(ACSC_cum3~FE_ACSC_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=ACSCsample)

c <- stargazer(mh, ph, acsc,
          apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),mdigits=2,
          add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(mh$fitted.values),4), 100*round(mean(ph$fitted.values),4), 100*round(mean(acsc$fitted.values),4))))

rbind(a[1:17], b[15:17], c[15:25])


#####################################################################
############             APPENDIX TABLE A2             ##############
#####################################################################

# Run all of the code in "Create_Analysis_Sample.R" with pcp_claims <- pcp_claims[PriorVisit==0] at the beginning
# Then run the code below

visit_composition <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;
SELECT
  PatientICN, COUNT(DISTINCT VisitDate) AS N_All
FROM(
SELECT DISTINCT
  cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate
FROM CDWWork.Outpat.Visit AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON visit.PrimaryStopCodeSID=sc.StopCodeSID
WHERE visit.VisitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND Visit.VisitDateTime<CAST('2018-03-01' AS datetime2(0)) AND StopCodeName NOT LIKE '%ADMIN%' AND StopCodeName NOT LIKE '%TELE%'
UNION(
SELECT
  cohort.PatientICN, CAST(AdmitDateTime AS DATE) AS VisitDate
FROM CDWWork.Inpat.Inpatient AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.AdmitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.AdmitDateTime)<=(365)
WHERE visit.AdmitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND visit.AdmitDateTime<CAST('2018-03-01' AS datetime2(0))
)
) AS z
GROUP BY PatientICN")

visit_composition <- data.table(visit_composition)
pcp_claims <- merge(x=pcp_claims, y=visit_composition, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_All[which(is.na(pcp_claims$N_All))] <- 0
pcp_claims$N_All[which(pcp_claims$N_All>=72)] <- 72

referrals <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT DISTINCT
  cohort.PatientICN,
  CASE WHEN con.ConsultSID IS NOT NULL THEN 1 ELSE 0 END AS Ref
FROM CDWWork.Con.Consult AS con
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON cohort.PatientICN=pat.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
INNER JOIN CDWWork.Dim.RequestService AS req
  ON con.ToRequestServiceSID=req.RequestServiceSID
LEFT OUTER JOIN CDWWork.Dim.AssociatedStopCode AS sc
  ON con.ToRequestServiceSID=sc.RequestServiceSID
INNER JOIN CDWWork.Dim.StopCode AS sc2
  ON sc.StopCodeSID=sc2.StopCodeSID
  WHERE RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND RequestDateTime<CAST('2018-06-08 00:00:00' AS datetime2(0))
  AND CPRSStatus NOT IN ('*%', 'DISCONTINUED', 'CANCELLED') 
  AND RequestType='C' -- a consult and not a procedure
  AND (AdministrativeFlag IS NULL OR NOT AdministrativeFlag='Y')
  AND InpatOutpat='O'
")

referrals <- data.table(referrals)
referrals <- referrals[PatientICN %in% pcp_claims$PatientICN]
referrals <- referrals[, .(Ref=max(Ref)), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=referrals, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$Ref[which(is.na(pcp_claims$Ref))] <- 0

pcp_claims$herc3 <- log(1+pcp_claims$herc_cum3)

hospital <- felm(PH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
death <- felm(Death3Year~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
Ref <- felm(Ref~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

mh <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
death <- felm(Death3Year~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
Ref <- felm(Ref~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

ph <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED or Hospitalizations", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
death <- felm(Death3Year~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
Ref <- felm(Ref~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

acsc <- stargazer(hospital, death, herc3, N_All, Ref,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, dep.var.labels=c("Circulatory ED or Hospitalizations", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

rbind(mh[1:17], ph[15:17], acsc[15:25])



########################################################################
##############             APPENDIX TABLE A3              ##############
########################################################################

ACSC <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; SELECT * FROM OMHO_QFR.ECON.appointments_acsc_nonVA")

ACSC$ACSCDate <- as.Date(ACSC$ACSCDate, "%Y-%m-%d")
ACSC <- merge(x=ACSC, y=pcp_claims[,c("PatientICN", "AppointmentRequestDate")], by.x="PatientICN", by.y="PatientICN")
ACSC$Post <- (ACSC$ACSCDate>ACSC$AppointmentRequestDate)

pcp_claims$ACSC_cum3_v3 <- (pcp_claims$PatientICN %in% ACSC$PatientICN[which(ACSC$Post==1)])
pcp_claims$ACSC_lag1_v3 <- (pcp_claims$PatientICN %in% ACSC$PatientICN[which(ACSC$Post==0)])

pcp_claims$ACSC_cum3_v3 <- pmax(pcp_claims$ACSC_cum3, pcp_claims$ACSC_cum3_v3)

pcp_claims$residual_MH <- resid(felm(MH_dummy_cum3_v3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))
pcp_claims$residual_PH <- resid(felm(PH_dummy_cum3_v3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))
pcp_claims$residual_ACSC <- resid(felm(ACSC_cum3_v3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))

drift <- pcp_claims[, `:=` (Sum_Resid_MH=sum(residual_MH, na.rm=T), 
                            Sum_Resid_PH=sum(residual_PH, na.rm=T),
                            Sum_Resid_ACSC=sum(residual_ACSC, na.rm=T),
                            AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]

drift$Z_MH_year <- (drift$Sum_Resid_MH - drift$residual_MH)/(drift$AnnualCaseCount - 1)
drift$Z_PH_year <- (drift$Sum_Resid_PH - drift$residual_PH)/(drift$AnnualCaseCount - 1)
drift$Z_ACSC_year <- (drift$Sum_Resid_ACSC - drift$residual_ACSC)/(drift$AnnualCaseCount - 1)

drift <- drift[,c("PatientICN", "MH_dummy_cum3_v3", "PH_dummy_cum3_v3", "ACSC_cum3_v3", "Z_MH_year", "Z_PH_year", "Z_ACSC_year", "StaffSSN", "AnnualCaseCount", "AssignYear")]

drift$N_bin <- "1"
drift$N_bin[which(drift$AnnualCaseCount>=10 & drift$AnnualCaseCount<25)] <- "2"
drift$N_bin[which(drift$AnnualCaseCount>=25 & drift$AnnualCaseCount<50)] <- "3"
drift$N_bin[which(drift$AnnualCaseCount>=50)] <- "4"

drift_wide <- reshape(drift, v.names=c("Z_MH_year", "Z_PH_year", "Z_ACSC_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide")
drift_wide <- data.frame(drift_wide)
drift_wide[is.na(drift_wide)] <- 0
drift_wide <- data.table(drift_wide)

drift <- merge(x=drift, y=drift_wide[,-c("PatientICN", "MH_dummy_cum3_v3", "PH_dummy_cum3_v3", "ACSC_cum3_v3")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

# Get 
drift$FE_MH_v3 <- lm(MH_dummy_cum3_v3~N_bin.2005*Z_MH_year.2005+N_bin.2006*Z_MH_year.2006+N_bin.2007*Z_MH_year.2007+N_bin.2008*Z_MH_year.2008+N_bin.2009*Z_MH_year.2009+N_bin.2010*Z_MH_year.2010+N_bin.2011*Z_MH_year.2011+N_bin.2012*Z_MH_year.2012+N_bin.2013*Z_MH_year.2013+N_bin.2014*Z_MH_year.2014+N_bin.2015*Z_MH_year.2015+N_bin.2016*Z_MH_year.2016+N_bin.2017*Z_MH_year.2017, data=drift)$fitted.values
drift$FE_PH_v3 <- lm(PH_dummy_cum3_v3~N_bin.2005*Z_PH_year.2005+N_bin.2006*Z_PH_year.2006+N_bin.2007*Z_PH_year.2007+N_bin.2008*Z_PH_year.2008+N_bin.2009*Z_PH_year.2009+N_bin.2010*Z_PH_year.2010+N_bin.2011*Z_PH_year.2011+N_bin.2012*Z_PH_year.2012+N_bin.2013*Z_PH_year.2013+N_bin.2014*Z_PH_year.2014+N_bin.2015*Z_PH_year.2015+N_bin.2016*Z_PH_year.2016+N_bin.2017*Z_PH_year.2017, data=drift)$fitted.values
drift$FE_ACSC_v3 <- lm(ACSC_cum3_v3~N_bin.2005*Z_ACSC_year.2005+N_bin.2006*Z_ACSC_year.2006+N_bin.2007*Z_ACSC_year.2007+N_bin.2008*Z_ACSC_year.2008+N_bin.2009*Z_ACSC_year.2009+N_bin.2010*Z_ACSC_year.2010+N_bin.2011*Z_ACSC_year.2011+N_bin.2012*Z_ACSC_year.2012+N_bin.2013*Z_ACSC_year.2013+N_bin.2014*Z_ACSC_year.2014+N_bin.2015*Z_ACSC_year.2015+N_bin.2016*Z_ACSC_year.2016+N_bin.2017*Z_ACSC_year.2017, data=drift)$fitted.values

pcp_claims <- merge(x=pcp_claims, y=drift[,c("PatientICN", "FE_MH_v3", "FE_PH_v3", "FE_ACSC_v3")], by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$FE_MH_std_v3 <- -scale(pcp_claims$FE_MH_v3)
pcp_claims$FE_PH_std_v3 <- -scale(pcp_claims$FE_PH_v3)
pcp_claims$FE_ACSC_std_v3 <- -scale(pcp_claims$FE_ACSC_v3)


hospital <- felm(PH_dummy_cum3_v1~FE_MH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death <- felm(Death3Year~FE_MH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_MH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_MH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref <- felm(Ref~FE_MH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

mh <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std_v3"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_PH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death <- felm(Death3Year~FE_PH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_PH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_PH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref <- felm(Ref~FE_PH_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

ph <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std_v3"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_ACSC_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death <- felm(Death3Year~FE_ACSC_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_ACSC_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_ACSC_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref <- felm(Ref~FE_ACSC_std_v3+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

acsc <- stargazer(hospital, death, herc3, N_All, Ref,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std_v3"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

rbind(mh[1:17], ph[15:17], acsc[15:25])



#####################################################################
############             APPENDIX TABLE A4             ##############
#####################################################################

# Run all of the code in "Create_Analysis_Sample.R" with pcp_claim <- pcp_claims[WaitTime>=0 & WaitTime<=14] at the beginning
# Then run the code below

visit_composition <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;
SELECT
  PatientICN, COUNT(DISTINCT VisitDate) AS N_All
FROM(
SELECT DISTINCT
  cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate
FROM CDWWork.Outpat.Visit AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON visit.PrimaryStopCodeSID=sc.StopCodeSID
WHERE visit.VisitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND Visit.VisitDateTime<CAST('2018-03-01' AS datetime2(0)) AND StopCodeName NOT LIKE '%ADMIN%' AND StopCodeName NOT LIKE '%TELE%'
UNION(
SELECT
  cohort.PatientICN, CAST(AdmitDateTime AS DATE) AS VisitDate
FROM CDWWork.Inpat.Inpatient AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.AdmitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.AdmitDateTime)<=(365)
WHERE visit.AdmitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND visit.AdmitDateTime<CAST('2018-03-01' AS datetime2(0))
)
) AS z
GROUP BY PatientICN")

visit_composition <- data.table(visit_composition)
pcp_claims <- merge(x=pcp_claims, y=visit_composition, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_All[which(is.na(pcp_claims$N_All))] <- 0
pcp_claims$N_All[which(pcp_claims$N_All>=72)] <- 72

referrals <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT DISTINCT
  cohort.PatientICN,
  CASE WHEN con.ConsultSID IS NOT NULL THEN 1 ELSE 0 END AS Ref
FROM CDWWork.Con.Consult AS con
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON cohort.PatientICN=pat.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
INNER JOIN CDWWork.Dim.RequestService AS req
  ON con.ToRequestServiceSID=req.RequestServiceSID
LEFT OUTER JOIN CDWWork.Dim.AssociatedStopCode AS sc
  ON con.ToRequestServiceSID=sc.RequestServiceSID
INNER JOIN CDWWork.Dim.StopCode AS sc2
  ON sc.StopCodeSID=sc2.StopCodeSID
  WHERE RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND RequestDateTime<CAST('2018-06-08 00:00:00' AS datetime2(0))
  AND CPRSStatus NOT IN ('*%', 'DISCONTINUED', 'CANCELLED') 
  AND RequestType='C' -- a consult and not a procedure
  AND (AdministrativeFlag IS NULL OR NOT AdministrativeFlag='Y')
  AND InpatOutpat='O'
")

referrals <- data.table(referrals)
referrals <- referrals[PatientICN %in% pcp_claims$PatientICN]
referrals <- referrals[, .(Ref=max(Ref)), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=referrals, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$Ref[which(is.na(pcp_claims$Ref))] <- 0

pcp_claims$herc3 <- log(1+pcp_claims$herc_cum3)

hospital <- felm(PH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
death <- felm(Death3Year~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
Ref <- felm(Ref~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

mh <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
death <- felm(Death3Year~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
Ref <- felm(Ref~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

ph <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED or Hospitalizations", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
death <- felm(Death3Year~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3 <- felm(herc3~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
Ref <- felm(Ref~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

acsc <- stargazer(hospital, death, herc3, N_All, Ref,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, dep.var.labels=c("Circulatory ED or Hospitalizations", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

rbind(mh[1:17], ph[15:17], acsc[15:25])


#########################################################
##########          APPENDIX TABLE A5          ##########
#########################################################

##### F-test for age for each facility
age_reg <- function(dt){
  model <- felm(Age~as.factor(StaffSSN)|YearXMonth+DayOfWeek+TimeOut_bins, data=dt)
  return(summary(model)$P.fstat[1])
}
pcp_claims[, `:=` (N_PCP=length(unique(StaffSSN))), by="InstitutionSID"]
ftest <- pcp_claims[N_PCP>2, .(p=age_reg(.SD)), by="InstitutionSID"]
random_facility <- pcp_claims[InstitutionSID %in% ftest$InstitutionSID[which(ftest$p>0.1)]]


hospital <- felm(PH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
death <- felm(Death3Year~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
herc3 <- felm(herc3~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
Ref <- felm(Ref~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)

mh <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
death <- felm(Death3Year~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
herc3 <- felm(herc3~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
Ref <- felm(Ref~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)

ph <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2)

hospital <- felm(PH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
death <- felm(Death3Year~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
herc3 <- felm(herc3~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility[DeathRelativeMonth>36])
N_All <- felm(N_All~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)
Ref <- felm(Ref~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=random_facility)

acsc <- stargazer(hospital, death, herc3, N_All, Ref,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(death$fitted.values),3), 100*round(mean(herc3$fitted.values),3), 100*round(mean(N_All$fitted.values),3), 100*round(mean(Ref$fitted.values),3))))

rbind(mh[1:17], ph[15:17], acsc[15:25])



############################################################################
################             APPENDIX TABLE A6             #################
############################################################################


visit_composition <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT
  PatientICN, COUNT(DISTINCT VisitDate) AS N_All
FROM(
SELECT DISTINCT
  cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate
FROM CDWWork.Outpat.Visit AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON visit.PrimaryStopCodeSID=sc.StopCodeSID
WHERE visit.VisitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND Visit.VisitDateTime<CAST('2018-03-01' AS datetime2(0)) AND StopCodeName NOT LIKE '%ADMIN%' AND StopCodeName NOT LIKE '%TELE%'
UNION(
SELECT
  cohort.PatientICN, CAST(AdmitDateTime AS DATE) AS VisitDate
FROM CDWWork.Inpat.Inpatient AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.AdmitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.AdmitDateTime)<=(365)
WHERE visit.AdmitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND visit.AdmitDateTime<CAST('2018-03-01' AS datetime2(0))
)
) AS z
GROUP BY PatientICN")

visit_composition <- data.table(visit_composition)
pcp_claims <- merge(x=pcp_claims, y=visit_composition, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_All[which(is.na(pcp_claims$N_All))] <- 0
pcp_claims$N_All[which(pcp_claims$N_All>=72)] <- 72

referrals <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT DISTINCT
  cohort.PatientICN,
  CASE WHEN con.ConsultSID IS NOT NULL THEN 1 ELSE 0 END AS Ref
FROM CDWWork.Con.Consult AS con
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON cohort.PatientICN=pat.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
INNER JOIN CDWWork.Dim.RequestService AS req
  ON con.ToRequestServiceSID=req.RequestServiceSID
LEFT OUTER JOIN CDWWork.Dim.AssociatedStopCode AS sc
  ON con.ToRequestServiceSID=sc.RequestServiceSID
INNER JOIN CDWWork.Dim.StopCode AS sc2
  ON sc.StopCodeSID=sc2.StopCodeSID
  WHERE RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND RequestDateTime<CAST('2018-06-08 00:00:00' AS datetime2(0))
  AND CPRSStatus NOT IN ('*%', 'DISCONTINUED', 'CANCELLED') 
  AND RequestType='C' -- a consult and not a procedure
  AND (AdministrativeFlag IS NULL OR NOT AdministrativeFlag='Y')
  AND InpatOutpat='O'
")

referrals <- data.table(referrals)
referrals <- referrals[PatientICN %in% pcp_claims$PatientICN]
referrals <- referrals[, .(Ref=max(Ref)), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=referrals, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$Ref[which(is.na(pcp_claims$Ref))] <- 0

pcp_claims$herc3 <- log(1+pcp_claims$herc_cum3)

hospital <- felm(PH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
death <- felm(Death3Year~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
herc3 <- felm(herc3~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)") & DeathRelativeMonth>36])
N_All <- felm(N_All~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
Ref <- felm(Ref~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])

mh <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

hospital <- felm(PH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
death <- felm(Death3Year~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
herc3 <- felm(herc3~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)") & DeathRelativeMonth>36])
N_All <- felm(N_All~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
Ref <- felm(Ref~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])

ph <- stargazer(hospital, death, herc3, N_All, Ref,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2)

hospital <- felm(PH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
death <- felm(Death3Year~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
herc3 <- felm(herc3~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)") & DeathRelativeMonth>36])
N_All <- felm(N_All~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])
Ref <- felm(Ref~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)")])

acsc <- stargazer(hospital, death, herc3, N_All, Ref,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(death$fitted.values),3), 100*round(mean(herc3$fitted.values),3), 100*round(mean(N_All$fitted.values),3), 100*round(mean(Ref$fitted.values),3))))

rbind(mh[1:17], ph[15:17], acsc[15:25])



#####################################################################
#############             APPENDIX TABLE A7             #############
#####################################################################

wrvu <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT DISTINCT 
  cohort.PatientICN, CAST(visit.VisitDateTime AS date) AS VisitDate, WRVU
FROM CDWWork.Outpat.Visit AS visit
LEFT OUTER JOIN(
  SELECT VisitSID, WRVU*Quantity AS WRVU
  FROM CDWWork.Outpat.VProcedure AS z
  INNER JOIN CDWWork.Dim.CPT AS cpt
    ON z.CPTSID=cpt.CPTSID
  INNER JOIN OMHO_QFR.ECON.wRVU_CPT AS wrvu
    ON cpt.CPTCode=wrvu.HCPCS
  WHERE VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND VisitDateTime<CAST('2018-06-01 00:00:00' AS datetime2(0))
) AS proce
  ON visit.VisitSID=proce.VisitSID
INNER JOIN CDWWork.Outpat.VProvider AS prov
  ON visit.VisitSID=prov.VisitSID
INNER JOIN CDWWork.SStaff.SStaff AS staff
  ON prov.ProviderSID=staff.StaffSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON staff.StaffSSN=cohort.StaffSSN AND pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime BETWEEN cohort.AppointmentRequestDate AND DATEADD(day, 365, cohort.AppointmentRequestDate)
WHERE visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND visit.VisitDateTime<CAST('2018-06-01 00:00:00' AS datetime2(0))
AND prov.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND prov.VisitDateTime<CAST('2018-06-01 00:00:00' AS datetime2(0))
")

wrvu <- data.table(wrvu)
wrvu <- wrvu[, .(WRVU=mean(WRVU)), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=wrvu, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$WRVU[which(is.na(pcp_claims$WRVU))] <- 0

MH <- felm(WRVU~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
           +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>12])
PH <- felm(WRVU~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
           +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>12])
ACSC <- felm(WRVU~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>12])

stargazer(MH, PH, ACSC,apply.coef=function(x){x*100}, apply.se=function(x){x*100},
          keep=c("FE_MH_std", "FE_PH_std", "FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
          digits=3)




#################################################################
############            APPENDIX TABLE A8            ############
#################################################################

### ED only
pcp_claims$residual_MH <- resid(felm(MH_dummy_ED_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))
pcp_claims$residual_PH <- resid(felm(PH_dummy_ED_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))

drift <- pcp_claims[, `:=` (Sum_Resid_MH=sum(residual_MH, na.rm=T), Sum_Resid_PH=sum(residual_PH, na.rm=T), AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]

drift$Z_MH_year <- (drift$Sum_Resid_MH - drift$residual_MH)/(drift$AnnualCaseCount - 1)
drift$Z_PH_year <- (drift$Sum_Resid_PH - drift$residual_PH)/(drift$AnnualCaseCount - 1)

drift <- drift[,c("PatientICN", "MH_dummy_ED_cum3", "PH_dummy_ED_cum3", "Z_MH_year", "Z_PH_year", "StaffSSN", "AnnualCaseCount", "AssignYear")]

drift$N_bin <- "1"
drift$N_bin[which(drift$AnnualCaseCount>=10 & drift$AnnualCaseCount<25)] <- "2"
drift$N_bin[which(drift$AnnualCaseCount>=25 & drift$AnnualCaseCount<50)] <- "3"
drift$N_bin[which(drift$AnnualCaseCount>=50)] <- "4"

drift_wide <- reshape(drift, v.names=c("Z_MH_year", "Z_PH_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide")
drift_wide <- data.frame(drift_wide)
drift_wide[is.na(drift_wide)] <- 0
drift_wide <- data.table(drift_wide)

drift <- merge(x=drift, y=drift_wide[,-c("PatientICN", "MH_dummy_ED_cum3", "PH_dummy_ED_cum3")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

drift$FE_MH_ED <- lm(MH_dummy_ED_cum3~N_bin.2005*Z_MH_year.2005+N_bin.2006*Z_MH_year.2006+N_bin.2007*Z_MH_year.2007+N_bin.2008*Z_MH_year.2008+N_bin.2009*Z_MH_year.2009+N_bin.2010*Z_MH_year.2010+N_bin.2011*Z_MH_year.2011+N_bin.2012*Z_MH_year.2012+N_bin.2013*Z_MH_year.2013+N_bin.2014*Z_MH_year.2014+N_bin.2015*Z_MH_year.2015+N_bin.2016*Z_MH_year.2016+N_bin.2017*Z_MH_year.2017, data=drift)$fitted.values
drift$FE_PH_ED <- lm(PH_dummy_ED_cum3~N_bin.2005*Z_PH_year.2005+N_bin.2006*Z_PH_year.2006+N_bin.2007*Z_PH_year.2007+N_bin.2008*Z_PH_year.2008+N_bin.2009*Z_PH_year.2009+N_bin.2010*Z_PH_year.2010+N_bin.2011*Z_PH_year.2011+N_bin.2012*Z_PH_year.2012+N_bin.2013*Z_PH_year.2013+N_bin.2014*Z_PH_year.2014+N_bin.2015*Z_PH_year.2015+N_bin.2016*Z_PH_year.2016+N_bin.2017*Z_PH_year.2017, data=drift)$fitted.values

pcp_claims <- merge(x=pcp_claims, y=drift[,c("PatientICN","FE_MH_ED", "FE_PH_ED")], by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$FE_MH_ED_std <- -scale(pcp_claims$FE_MH_ED)
pcp_claims$FE_PH_ED_std <- -scale(pcp_claims$FE_PH_ED)

# Create Value Add
ValueAdd <- pcp_claims[, .(FirstYear=min(as.numeric(VisitYear)), LastYear=max(as.numeric(VisitYear)), FrequentYear=as.numeric(names(sort(table(VisitYear), dec=T)[1])), N=length(VisitSID),
                           ProviderType=names(sort(table(ProviderType), dec=T)[1]), Sta3n=names(sort(table(Sta3n), dec=T)[1]), InstitutionSID=names(sort(table(InstitutionSID), dec=T)[1])), by="StaffSSN"]
ValueAdd$StaffSSN <- as.numeric(as.character(ValueAdd$StaffSSN))

FE <- drift[, .(FE_MH_ED=mean(FE_MH_ED), FE_PH_ED=mean(FE_PH_ED), AnnualCaseCount=length(PatientICN)), by="StaffSSN"]
FE$StaffSSN <- as.numeric(as.character(FE$StaffSSN))
ValueAdd <- merge(x=ValueAdd, y=FE, by.x="StaffSSN", by.y="StaffSSN")

### Inpatient only
pcp_claims$residual_MH <- resid(felm(MH_dummy_HOSP_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))
pcp_claims$residual_PH <- resid(felm(PH_dummy_HOSP_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))

drift <- pcp_claims[, `:=` (Sum_Resid_MH=sum(residual_MH, na.rm=T), Sum_Resid_PH=sum(residual_PH, na.rm=T), AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]

drift$Z_MH_year <- (drift$Sum_Resid_MH - drift$residual_MH)/(drift$AnnualCaseCount - 1)
drift$Z_PH_year <- (drift$Sum_Resid_PH - drift$residual_PH)/(drift$AnnualCaseCount - 1)

drift <- drift[,c("PatientICN", "MH_dummy_HOSP_cum3", "PH_dummy_HOSP_cum3", "Z_MH_year", "Z_PH_year", "StaffSSN", "AnnualCaseCount", "AssignYear")]

drift$N_bin <- "1"
drift$N_bin[which(drift$AnnualCaseCount>=10 & drift$AnnualCaseCount<25)] <- "2"
drift$N_bin[which(drift$AnnualCaseCount>=25 & drift$AnnualCaseCount<50)] <- "3"
drift$N_bin[which(drift$AnnualCaseCount>=50)] <- "4"

drift_wide <- reshape(drift, v.names=c("Z_MH_year", "Z_PH_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide")
drift_wide <- data.frame(drift_wide)
drift_wide[is.na(drift_wide)] <- 0
drift_wide <- data.table(drift_wide)

drift <- merge(x=drift, y=drift_wide[,-c("PatientICN", "MH_dummy_HOSP_cum3", "PH_dummy_HOSP_cum3")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

drift$FE_MH_HOSP <- lm(MH_dummy_HOSP_cum3~N_bin.2005*Z_MH_year.2005+N_bin.2006*Z_MH_year.2006+N_bin.2007*Z_MH_year.2007+N_bin.2008*Z_MH_year.2008+N_bin.2009*Z_MH_year.2009+N_bin.2010*Z_MH_year.2010+N_bin.2011*Z_MH_year.2011+N_bin.2012*Z_MH_year.2012+N_bin.2013*Z_MH_year.2013+N_bin.2014*Z_MH_year.2014+N_bin.2015*Z_MH_year.2015+N_bin.2016*Z_MH_year.2016+N_bin.2017*Z_MH_year.2017, data=drift)$fitted.values
drift$FE_PH_HOSP <- lm(PH_dummy_HOSP_cum3~N_bin.2005*Z_PH_year.2005+N_bin.2006*Z_PH_year.2006+N_bin.2007*Z_PH_year.2007+N_bin.2008*Z_PH_year.2008+N_bin.2009*Z_PH_year.2009+N_bin.2010*Z_PH_year.2010+N_bin.2011*Z_PH_year.2011+N_bin.2012*Z_PH_year.2012+N_bin.2013*Z_PH_year.2013+N_bin.2014*Z_PH_year.2014+N_bin.2015*Z_PH_year.2015+N_bin.2016*Z_PH_year.2016+N_bin.2017*Z_PH_year.2017, data=drift)$fitted.values

pcp_claims <- merge(x=pcp_claims, y=drift[,c("PatientICN","FE_MH_HOSP", "FE_PH_HOSP")], by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$FE_MH_HOSP_std <- -scale(pcp_claims$FE_MH_HOSP)
pcp_claims$FE_PH_HOSP_std <- -scale(pcp_claims$FE_PH_HOSP)


FE <- drift[, .(FE_MH_HOSP=mean(FE_MH_HOSP), FE_PH_HOSP=mean(FE_PH_HOSP), AnnualCaseCount=length(PatientICN)), by="StaffSSN"]
FE$StaffSSN <- as.numeric(as.character(FE$StaffSSN))
ValueAdd <- merge(x=ValueAdd, y=FE, by.x="StaffSSN", by.y="StaffSSN")

# Standardize and take negative
ValueAdd$FE_MH_ED_std <- -scale(ValueAdd$FE_MH_ED)
ValueAdd$FE_PH_ED_std <- -scale(ValueAdd$FE_PH_ED)
ValueAdd$FE_MH_HOSP_std <- -scale(ValueAdd$FE_MH_HOSP)
ValueAdd$FE_PH_HOSP_std <- -scale(ValueAdd$FE_PH_HOSP)

# What is correlation between ED and Hospital effectiveness?
summary(felm(FE_MH_ED_std~0+FE_MH_HOSP_std|Sta3n, data=ValueAdd, weights=ValueAdd$N)); summary(felm(FE_PH_ED_std~0+FE_PH_HOSP_std|Sta3n, data=ValueAdd, weights=ValueAdd$N))


visit_composition <- sqlQuery(db.handle,  "set nocount on; set ansi_warnings off;
SELECT
  PatientICN, COUNT(DISTINCT VisitDate) AS N_All
FROM(
SELECT DISTINCT
  cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate
FROM CDWWork.Outpat.Visit AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON visit.PrimaryStopCodeSID=sc.StopCodeSID
WHERE visit.VisitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND Visit.VisitDateTime<CAST('2018-03-01' AS datetime2(0)) AND StopCodeName NOT LIKE '%ADMIN%' AND StopCodeName NOT LIKE '%TELE%'
UNION(
SELECT
  cohort.PatientICN, CAST(AdmitDateTime AS DATE) AS VisitDate
FROM CDWWork.Inpat.Inpatient AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.AdmitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.AdmitDateTime)<=(365)
WHERE visit.AdmitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND visit.AdmitDateTime<CAST('2018-03-01' AS datetime2(0))
)
) AS z
GROUP BY PatientICN")

visit_composition <- data.table(visit_composition)
pcp_claims <- merge(x=pcp_claims, y=visit_composition, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_All[which(is.na(pcp_claims$N_All))] <- 0
pcp_claims$N_All[which(pcp_claims$N_All>=72)] <- 72

referrals <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;
SELECT DISTINCT
  cohort.PatientICN,
  CASE WHEN con.ConsultSID IS NOT NULL THEN 1 ELSE 0 END AS Ref
FROM CDWWork.Con.Consult AS con
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON cohort.PatientICN=pat.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
INNER JOIN CDWWork.Dim.RequestService AS req
  ON con.ToRequestServiceSID=req.RequestServiceSID
LEFT OUTER JOIN CDWWork.Dim.AssociatedStopCode AS sc
  ON con.ToRequestServiceSID=sc.RequestServiceSID
INNER JOIN CDWWork.Dim.StopCode AS sc2
  ON sc.StopCodeSID=sc2.StopCodeSID
  WHERE RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND RequestDateTime<CAST('2018-06-08 00:00:00' AS datetime2(0))
  AND CPRSStatus NOT IN ('*%', 'DISCONTINUED', 'CANCELLED') 
  AND RequestType='C' -- a consult and not a procedure
  AND (AdministrativeFlag IS NULL OR NOT AdministrativeFlag='Y')
  AND InpatOutpat='O'
")

referrals <- data.table(referrals)
referrals <- referrals[PatientICN %in% pcp_claims$PatientICN]
referrals <- referrals[, .(Ref=max(Ref)), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=referrals, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$Ref[which(is.na(pcp_claims$Ref))] <- 0

pcp_claims$herc3 <- log(1+pcp_claims$herc_cum3)

# MH
hospital_ED <- felm(PH_dummy_cum3_v1~FE_MH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                    +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hospital_Hosp <- felm(PH_dummy_cum3_v1~FE_MH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                      +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death_ED <- felm(Death3Year~FE_MH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death_Hosp <- felm(Death3Year~FE_MH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3_ED <- felm(herc3~FE_MH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
herc3_HOSP <- felm(herc3~FE_MH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All_ED <- felm(N_All~FE_MH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_All_HOSP <- felm(N_All~FE_MH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_ED <- felm(Ref~FE_MH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_HOSP <- felm(Ref~FE_MH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

mh <- stargazer(hospital_ED, hospital_Hosp, death_ED, death_Hosp, herc3_ED, herc3_HOSP, N_All_ED, N_All_HOSP, Ref_ED, Ref_HOSP,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_ED_std", "FE_MH_HOSP_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

# PH
hospital_ED <- felm(PH_dummy_cum3_v1~FE_PH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                    +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hospital_Hosp <- felm(PH_dummy_cum3_v1~FE_PH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                      +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death_ED <- felm(Death3Year~FE_PH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
death_Hosp <- felm(Death3Year~FE_PH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc3_ED <- felm(herc3~FE_PH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
herc3_HOSP <- felm(herc3~FE_PH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])
N_All_ED <- felm(N_All~FE_PH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_All_HOSP <- felm(N_All~FE_PH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_ED <- felm(Ref~FE_PH_ED_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_HOSP <- felm(Ref~FE_PH_HOSP_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                 +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

ph <- stargazer(hospital_ED, hospital_Hosp, death_ED, death_Hosp, herc3_ED, herc3_HOSP, N_All_ED, N_All_HOSP, Ref_ED, Ref_HOSP,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_ED_std", "FE_PH_HOSP_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Circulatory ED/Hosp", "3 Year Mortality", "Log 3Y Avg Cost", "Number of Visits", "Any Referrals"))

rbind(mh[1:17], ph[15:25])



########################################################################
################          APPENDIX TABLE A10            ################
########################################################################

# Providers with at least 20 cases each period
pcp_claims <- pcp_claims[, `:=` (N_0511=length(VisitSID[which(AssignYear>=2005 & AssignYear<=2011)]), N_1217=length(VisitSID[which(AssignYear>=2012 & AssignYear<=2017)])), by="StaffSSN"]
pcp_claims <- pcp_claims[N_0511>=20 & N_1217>=20] # Drops us from 7550 to 2566 providers

sample_0511 <- pcp_claims[AssignYear>=2005 & AssignYear<=2011]
sample_1217 <- pcp_claims[AssignYear>=2012 & AssignYear<=2017]

### 2005-2011
sample_0511$residual_MH <- resid(felm(MH_dummy_cum3_v1~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=sample_0511))
sample_0511$residual_PH <- resid(felm(PH_dummy_cum3_v1~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=sample_0511))
sample_0511$residual_ACSC <- resid(felm(ACSC_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=sample_0511))

drift <- sample_0511[, `:=` (Sum_Resid_MH=sum(residual_MH, na.rm=T), 
                             Sum_Resid_PH=sum(residual_PH, na.rm=T),
                             Sum_Resid_ACSC=sum(residual_ACSC, na.rm=T),
                             AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]

drift$Z_MH_year <- (drift$Sum_Resid_MH - drift$residual_MH)/(drift$AnnualCaseCount - 1)
drift$Z_PH_year <- (drift$Sum_Resid_PH - drift$residual_PH)/(drift$AnnualCaseCount - 1)
drift$Z_ACSC_year <- (drift$Sum_Resid_ACSC - drift$residual_ACSC)/(drift$AnnualCaseCount - 1)

drift <- drift[,c("PatientICN", "MH_dummy_cum3_v1", "PH_dummy_cum3_v1", "ACSC_cum3", "Z_MH_year", "Z_PH_year", "Z_ACSC_year", "StaffSSN", "AnnualCaseCount", "AssignYear")]

drift$N_bin <- "1"
drift$N_bin[which(drift$AnnualCaseCount>=10 & drift$AnnualCaseCount<25)] <- "2"
drift$N_bin[which(drift$AnnualCaseCount>=25 & drift$AnnualCaseCount<50)] <- "3"
drift$N_bin[which(drift$AnnualCaseCount>=50)] <- "4"

drift_wide <- reshape(drift, v.names=c("Z_MH_year", "Z_PH_year", "Z_ACSC_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide")
drift_wide <- data.frame(drift_wide)
drift_wide[is.na(drift_wide)] <- 0
drift_wide <- data.table(drift_wide)

drift <- merge(x=drift, y=drift_wide[,-c("PatientICN", "MH_dummy_cum3_v1", "PH_dummy_cum3_v1", "ACSC_cum3")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

drift$FE_MH <- lm(MH_dummy_cum3_v1~N_bin.2005*Z_MH_year.2005+N_bin.2006*Z_MH_year.2006+N_bin.2007*Z_MH_year.2007+N_bin.2008*Z_MH_year.2008+N_bin.2009*Z_MH_year.2009+N_bin.2010*Z_MH_year.2010+N_bin.2011*Z_MH_year.2011, data=drift)$fitted.values
drift$FE_PH <- lm(PH_dummy_cum3_v1~N_bin.2005*Z_PH_year.2005+N_bin.2006*Z_PH_year.2006+N_bin.2007*Z_PH_year.2007+N_bin.2008*Z_PH_year.2008+N_bin.2009*Z_PH_year.2009+N_bin.2010*Z_PH_year.2010+N_bin.2011*Z_PH_year.2011, data=drift)$fitted.values
drift$FE_ACSC <- lm(ACSC_cum3~N_bin.2005*Z_ACSC_year.2005+N_bin.2006*Z_ACSC_year.2006+N_bin.2007*Z_ACSC_year.2007+N_bin.2008*Z_ACSC_year.2008+N_bin.2009*Z_ACSC_year.2009+N_bin.2010*Z_ACSC_year.2010+N_bin.2011*Z_ACSC_year.2011, data=drift)$fitted.values

### Create Value Add
FE_0511 <- drift[, .(FE_MH_0511=-mean(FE_MH), FE_PH_0511=-mean(FE_PH), FE_ACSC_0511=-mean(FE_ACSC)), by="StaffSSN"]
FE_0511$StaffSSN <- as.numeric(as.character(FE_0511$StaffSSN))


### 2012-2017
sample_1217$residual_MH <- resid(felm(MH_dummy_cum3_v1~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=sample_1217))
sample_1217$residual_PH <- resid(felm(PH_dummy_cum3_v1~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=sample_1217))
sample_1217$residual_ACSC <- resid(felm(ACSC_cum3~Black+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=sample_1217))

drift <- sample_1217[, `:=` (Sum_Resid_MH=sum(residual_MH, na.rm=T), 
                             Sum_Resid_PH=sum(residual_PH, na.rm=T),
                             Sum_Resid_ACSC=sum(residual_ACSC, na.rm=T),
                             AnnualCaseCount=length(PatientICN)), by=c("StaffSSN", "AssignYear")]

drift$Z_MH_year <- (drift$Sum_Resid_MH - drift$residual_MH)/(drift$AnnualCaseCount - 1)
drift$Z_PH_year <- (drift$Sum_Resid_PH - drift$residual_PH)/(drift$AnnualCaseCount - 1)
drift$Z_ACSC_year <- (drift$Sum_Resid_ACSC - drift$residual_ACSC)/(drift$AnnualCaseCount - 1)

drift <- drift[,c("PatientICN", "MH_dummy_cum3_v1", "PH_dummy_cum3_v1", "ACSC_cum3", "Z_MH_year", "Z_PH_year", "Z_ACSC_year", "StaffSSN", "AnnualCaseCount", "AssignYear")]

drift$N_bin <- "1"
drift$N_bin[which(drift$AnnualCaseCount>=10 & drift$AnnualCaseCount<25)] <- "2"
drift$N_bin[which(drift$AnnualCaseCount>=25 & drift$AnnualCaseCount<50)] <- "3"
drift$N_bin[which(drift$AnnualCaseCount>=50)] <- "4"

drift_wide <- reshape(drift, v.names=c("Z_MH_year", "Z_PH_year", "Z_ACSC_year", "N_bin"), idvar="StaffSSN", timevar="AssignYear", direct="wide")
drift_wide <- data.frame(drift_wide)
drift_wide[is.na(drift_wide)] <- 0
drift_wide <- data.table(drift_wide)

drift <- merge(x=drift, y=drift_wide[,-c("PatientICN", "MH_dummy_cum3_v1", "PH_dummy_cum3_v1", "ACSC_cum3")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

drift$FE_MH <- lm(MH_dummy_cum3_v1~N_bin.2012*Z_MH_year.2012+N_bin.2013*Z_MH_year.2013+N_bin.2014*Z_MH_year.2014+N_bin.2015*Z_MH_year.2015+N_bin.2016*Z_MH_year.2016+N_bin.2017*Z_MH_year.2017, data=drift)$fitted.values
drift$FE_PH <- lm(PH_dummy_cum3_v1~N_bin.2012*Z_PH_year.2012+N_bin.2013*Z_PH_year.2013+N_bin.2014*Z_PH_year.2014+N_bin.2015*Z_PH_year.2015+N_bin.2016*Z_PH_year.2016+N_bin.2017*Z_PH_year.2017, data=drift)$fitted.values
drift$FE_ACSC <- lm(ACSC_cum3~N_bin.2012*Z_ACSC_year.2012+N_bin.2013*Z_ACSC_year.2013+N_bin.2014*Z_ACSC_year.2014+N_bin.2015*Z_ACSC_year.2015+N_bin.2016*Z_ACSC_year.2016+N_bin.2017*Z_ACSC_year.2017, data=drift)$fitted.values

### Create Value Add
FE_1217 <- drift[, .(FE_MH_1217=-mean(FE_MH), FE_PH_1217=-mean(FE_PH), FE_ACSC_1217=-mean(FE_ACSC)), by="StaffSSN"]
FE_1217$StaffSSN <- as.numeric(as.character(FE_1217$StaffSSN))


FE <- merge(x=FE_0511, y=FE_1217, by.x="StaffSSN", by.y="StaffSSN")
FE$MH_pre_rank <- order(FE$FE_MH_0511, decreasing=TRUE)
FE$MH_post_rank <- order(FE$FE_MH_1217, decreasing=TRUE)
FE$PH_pre_rank <- order(FE$FE_PH_0511, decreasing=TRUE)
FE$PH_post_rank <- order(FE$FE_PH_1217, decreasing=TRUE)
FE$ACSC_pre_rank <- order(FE$FE_ACSC_0511, decreasing=TRUE)
FE$ACSC_post_rank <- order(FE$FE_ACSC_1217, decreasing=TRUE)

mh <- lm(FE_MH_1217~0+FE_MH_0511, data=FE)
ph <- lm(FE_PH_1217~0+FE_PH_0511, data=FE)
acsc <- lm(FE_ACSC_1217~0+FE_ACSC_0511, data=FE)

stargazer(mh, ph, acsc, keep=c("FE_MH_0511", "FE_PH_0511", "FE_ACSC_0511"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"), digits=2)


##############################################################
#############          APPENDIX TABLE A11        #############
##############################################################

# MH
MH <- round(c(
  o_delta(y="MH_dummy_cum3_v1", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="MH_dummy_cum3_v1", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="PH_dummy_cum3_v1", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="PH_dummy_cum3_v1", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="ACSC_cum3", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="ACSC_cum3", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="Death3Year", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="Death3Year", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="herc3", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="herc3", x="FE_MH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1]
),2)

# Circulatory
circulatory <- round(c(
  o_delta(y="MH_dummy_cum3_v1", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="MH_dummy_cum3_v1", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="PH_dummy_cum3_v1", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="PH_dummy_cum3_v1", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="ACSC_cum3", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="ACSC_cum3", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="Death3Year", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="Death3Year", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="herc3", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="herc3", x="FE_PH_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1]
), 2)

# ACSC
ACSC <- round(c(
  o_delta(y="MH_dummy_cum3_v1", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="MH_dummy_cum3_v1", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="PH_dummy_cum3_v1", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="PH_dummy_cum3_v1", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="ACSC_cum3", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="ACSC_cum3", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="Death3Year", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="Death3Year", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1],
  o_delta(y="herc3", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.15)$Value[1],
  o_delta(y="herc3", x="FE_ACSC_std", con="RaceEthnicity+Married+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime+EnrollmentPriority", m="InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins+Age_bins", data=pcp_claims, type="lm",R2max=0.20)$Value[1]
), 2)

rbind(MH, circulatory, ACSC)
